###############################################
#Jordan Mancl 12-01-2022                      #
#Measure D1-D4, D2-D3 distances               #
#Measure D1-D2-D3-D4 dihedral angle           #
#Run from VMD Tkconsole                       #
#Output single file for easy transfer to Excel#
###############################################

#define domains
set ad1 [atomselect top "chain A and resid 51 to 285"]
set ad2 [atomselect top "chain A and resid 286 to 515"]
set ad3 [atomselect top "chain A and resid 542 to 768"]
set ad4 [atomselect top "chain A and (resid 769 to 963 or resid 989 to 1007)"]
set bd1 [atomselect top "chain B and resid 51 to 285"]
set bd2 [atomselect top "chain B and resid 286 to 515"]
set bd3 [atomselect top "chain B and resid 542 to 768"]
set bd4 [atomselect top "chain B and (resid 769 to 963 or resid 989 to 1007)"]

#get frames in simulation
set nf [molinfo top get numframes]

#open output file
set outfile [open measurements.dat w]

#write column headings
puts $outfile "Frame AD1-D4 AD2-D3 BD1-D4 BD2-D3 DihedralA DihedralB"

#loop over all frames
for {set i 0} {$i < $nf} {incr i} {

	#write out frame number and update selections to current frame
	puts "frame $i of $nf"
	$ad1 frame $i
	$ad2 frame $i
	$ad3 frame $i
	$ad4 frame $i
	$bd1 frame $i
	$bd2 frame $i
	$bd3 frame $i
	$bd4 frame $i

	#find COM for each selection
	set com_ad1 [measure center $ad1 weight mass]
	puts "com_ad1 $com_ad1"
	set com_ad2 [measure center $ad2 weight mass]
	puts "com_ad2 $com_ad2"
	set com_ad3 [measure center $ad3 weight mass]
	puts "com_ad3 $com_ad3"
	set com_ad4 [measure center $ad4 weight mass]
	puts "com_ad4 $com_ad4"
	set com_bd1 [measure center $bd1 weight mass]
	puts "com_bd1 $com_bd1"
	set com_bd2 [measure center $bd2 weight mass]
	puts "com_bd2 $com_bd2"
	set com_bd3 [measure center $bd3 weight mass]
	puts "com_bd3 $com_bd3"
	set com_bd4 [measure center $bd4 weight mass]
	puts "com_bd4 $com_bd4"

	#at each frame i find the chain A d1-d4 distance, assign value to array
	set ad1d4($i.r) [veclength [vecsub $com_ad1 $com_ad4]]

	#now do chain A d2-d3 distance
	set ad2d3($i.r) [veclength [vecsub $com_ad2 $com_ad3]]

	#now do chain B d1-d4 distance
	set bd1d4($i.r) [veclength [vecsub $com_bd1 $com_bd4]]

	#now do chain B d2-d3 distance
	set bd2d3($i.r) [veclength [vecsub $com_bd2 $com_bd3]]

	## CALCULATE DIHEDRAL ANGLES
	# n1 and n2 = normal vectors formed by atoms 1 2 3 and 2 3 4
	# n1 = (x1 y1 z1)
	# n2 = (x2 y2 z2)
	# angle = acos (n1 . n2)/(|n1|*|n2|)
	#
	# Given 3 points: (x1 y1 z1) (x2 y2 z2) (x3 y3 z3)
	# Two vectors are: (x2-x1 y2-y1 z2-z1) (x3-x2 y3-y2 z3-z2)
	# Normal Vector = | i j k |
	# | x1-x2 y1-y2 z1-z2 |
	# | x2-x3 y2-y3 z2-z3 |
	#

	#chain A dihedral
	set a12 [vecsub $com_ad2 $com_ad1]
	set a21 [vecsub $com_ad1 $com_ad2]
	set a23 [vecsub $com_ad3 $com_ad2]
	set a32 [vecsub $com_ad2 $com_ad3]
	set a34 [vecsub $com_ad4 $com_ad3]
	set n1 [veccross $a12 $a23]
	set n2 [veccross $a23 $a34]
	set cosangle_a [expr [vecdot $n1 $n2]/([veclength $n1]*[veclength $n2])]
	set angle_a($i.r) [expr {180 * acos($cosangle_a) / 3.14159265358}]
	set sign_a [vecdot [veccross $n1 $n2] $a23]
	if {$sign_a < 0 } {
		set angle_a($i.r) [expr {0 - $angle_a($i.r)}]
	} elseif {$sign_a == 0 } {
		set angle_a($i.r) 180.0
	}
	puts $angle_a($i.r)

	#chain B dihedral
                set b12 [vecsub $com_bd2 $com_bd1]
                set b21 [vecsub $com_bd1 $com_bd2]
                set b23 [vecsub $com_bd3 $com_bd2]
                set b32 [vecsub $com_bd2 $com_bd3]
                set b34 [vecsub $com_bd4 $com_bd3]
                set n3 [veccross $b12 $b23]
                set n4 [veccross $b23 $b34]
                set cosangle_b [expr [vecdot $n3 $n4]/([veclength $n3]*[veclength $n4])]
                set angle_b($i.r) [expr {180 * acos($cosangle_b) / 3.14159265358}]
                set sign_b [vecdot [veccross $n3 $n4] $b23]
                if {$sign_b < 0 } {
                        set angle_b($i.r) [expr {0 - $angle_b($i.r)}]
                } elseif {$sign_b == 0 } {
                        set angle_b($i.r) 180.0
                }
                puts $angle_b($i.r)


	#write values to output file
	puts $outfile "$i $ad1d4($i.r) $ad2d3($i.r) $bd1d4($i.r) $bd2d3($i.r) $angle_a($i.r) $angle_b($i.r)"

	#end loop
}

#close the output file
close $outfile
